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Abstract 

interpolation  scheme  is  described  for  values  of  position,  gradient  and  Hessian  at 
scattered  points  in  three  variables.  The  domain  is  assumed  to  have  been  tesselated  into 
tetrahedra.  The  interpolant  has  local  support,  is  globally  once  differentiable,  piecewise 
polynomial,  and  reproduces  polynomials  of  degree  up  to  three  exactly.  The  scheme  has 
been  implemented  in  a  FORTRAN  research  code. 
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SIGNIFICANCE  AND  EXPLANATION 


This  report  describes  an  interpolation  scheme  for  scattered  data  in  three  independent 
variables,  as  they  may  arise,  e.g.  in  the  desciption  of  temperature  or  pressure  in  a  three- 
dimensional  solid.  The  interpolant  will  be  once  differentiable.  The  user  must  supply 
function  values  at  the  scattered  points  and  also  values  of  first  and  second  order  derivatives. 
The  domain  of  the  interpolant  is  the  convex  hull  of  the  given  points  and  is  assumed  to 
have  been  tesselated  with  tetrahedra.  The  scheme  is  piecewise  polynomial,  and  reproduces 
all  polynomials  of  degree  up  to  three  exactly.  A  non-portable  FORTRAN  research  code 
is  available  from  the  author. 
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Appendix:  Coefficients  of  the  Interpolant 
1.  Introduction 

Multivariate  interpolation  problems  arise  in  many  situations  where  physical  phenom¬ 
ena  involving  several  space  variables  (and  possibly  time  as  well)  are  modeled,  designed, 
analyzed,  or  simulated.  One  application  might  be  temperature  as  a  (trivariate)  function 
of  the  three  spatial  variables.  The  analyst  is  given  a  set  of  four-dimensional  data  through 
which  he  wishes  to  pass  a  smooth  trivariate  function.  In  a  design  context,  he  may  wish 
to  alter  the  data  interactively  until  a  satisfactory  function  has  been  found.  Because  the 
geometric  information  may  be  located  arbitrarily  in  four-dimensional  space,  the  interpo¬ 
lation  scheme  must  be  able  to  handle  such  data.  A  frequently  used  approach  consists  of 
smoothly  piecing  together  interpolants  defined  only  on  parts  of  the  underlying  domain. 
Such  schemes  have  the  property  of  being  local,  i.e.  the  value  of  the  interpolant  at  a  given 
point  is  dependent  only  upon  a  small  amount  of  data  close  to  the  point  of  evaluation. 
This  has  two  important  practical  advantages:  Firstly,  the  cost  of  evaluation  of  a  scheme  is 
largely  independent  of  the  amount  of  data,  and,  secondly,  a  change  in  the  data  alters  the 
interpolant  only  in  a  small  region  around  those  points  where  the  change  takes  place. 

We  will  assume  that  the  domain  of  interest  has  been  tesselated  into  tetrahedra,  with 
the  vertices  of  the  tetrahedra  being  the  data  points.  How  the  tesselation  may  be  accom¬ 
plished  is  an  important  and  difficult  question  which  which  will  not  be  addressed  in  this 
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paper.  However,  see  Barnhill  and  Little,  1984,  for  some  answers  to  this  question.  An 
interpolant  that  is  defined  piecewise  on  tetrahedra  is  called  tetrahedral. 

Important  attributes  of  a  scheme  include  the  maximum  degree  of  the  derivatives  that 
have  to  be  supplied  at  the  data  points  and  the  degree  of  differentiability  of  the  interpolant. 
We  use  the  term  C  scheme  for  C®  data  to  describe  a  globally  r-times  differentiable  scheme 
that  requires  derivatives  through  s-th  order  as  data.  Usually,  a  user  requires  his  scheme 
to  be  at  least  C',  but  is  unable  to  provide  data  of  degree  larger  than  0. 

Several  tetrahedral  schemes  are  known.  Alfeld,  1984a,  describes  a  C'  scheme  for  C’ 
data.  Barnhill  and  Little,  1984,  describe  a  similar  scheme  that  requires  some  auxilliary 
data.  In  another  paper,  Alfeld,  1984b,  describes  a  class  of  C""  schemes  for  C'"  data,  for 
arbitrary  m  and  in  arbitrarily  many  variables. 

However,  all  of  the  existing  tetrahedral  schemes  yield  rational  interpolants.  One  would 
like  to  have  polynomial  interpolants  for  several  reasons.  Firstly,  they  are  more  efficiently 
evaluated.  Secondly,  they  are  more  easily  integrated.  This  is  useful,  for  example,  in 
the  finite  element  technique  for  the  solution  of  partial  differential  equations  or  in  certain 
techniques  for  the  generation  of  missing  derivative  data  (see  Alfeld, 1984c). 

In  this  paper,  we  derive  a  C’  scheme  for  data  which  yields  an  interpolant  that 
is  piecewise  quintic  on  each  tetrahedron.  The  scheme  is  modeled  after  the  well-known 
bivariate  Clough-Tocher  scheme  (Strang  and  Fix,  1973,  p.  82).  That  scheme  is  a  C* 
scheme  for  C*  data.  Our  scheme  requires  data.  Attempts  to  obtain  a  trivariate 
Clough-Tocher  type  scheme  for  C*  data  have  been  undertaken  by  the  author  and  others, 
but  to  date  have  been  unsuccessful. 

However,  even  requiring  data  is  a  significant  gain  over  what  would  be  required  for 
straight  polynomial  interpolation  on  a  tetrahedron.  !^eni§ek,  1973,  has  shown  that  such  a 
scheme  would  require  a  polynomial  of  degree  at  least  9  (with  220  parameters)  and  at  least 
C*  data. 

The  paper  is  organized  as  follows:  In  section  2,  an  explicit  expression  for  the  bivariate 
Clough-Tocher  scheme  is  derived.  The  purpose  of  that  section  is  to  introduce  the  appro¬ 
priate  machinery  in  a  familiar  context.  Although  the  bivariate  scheme  has  been  used  for 
a  long  time,  to  the  author’s  knowledge  this  is  the  first  time  that  explicit  expressions  are 
given.  However,  a  reader  familiar  with  barycentric  coordinates  may  wish  to  read  only 
step  0  in  section  2  and  then  skip  to  section  three  which  introduces  the  trivariate  scheme. 
Section  4  briefly  discusses  computational  aspects. 

2.  The  Bivariate  Clough-Tocher  Scheme 

The  derivation  of  the  trivariate  Clough-Tocher  scheme  is  quite  involved  algebraically, 
and  difficult  to  illustrate  graphically  because  of  the  three-dimensional  domain.  Therefore, 
in  this  section,  the  well-known  bivariate  Clough-Tocher  scheme  is  derived  using  precisely 
the  same  approach  that  will  also  be  employed  for  the  derivation  of  the  trivariate  scheme. 
The  next  section  will  then  amount  to  a  description  of  the  differences  between  the  bivariate 
and  trivariate  schemes. 

The  bivariate  Clough-Tocher  scheme  uses  as  data  position  and  gradient  at  the  vertices 
of  a  general  triangle.  The  triangle  is  divided  about  its  centroid  into  three  subtriangles. 
We  refer  to  the  overall  triangle  as  the  macrotriangle  and  to  its  subtriangles  as  microtri- 
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angles.  The  interpolant  is  constructed  to  be  a  cubic  polynomial  on  each  microtriangle. 
It  is  differentiable  over  the  macrotriangle.  Continuity  between  macrotriangles  is  ensured 
by  interpolation  to  the  data.  Differentiability  between  macrotriangles  is  forced  by  the 
requirement  that  the  first  order  perpendicular  cross-boundary  derivative  along  edges  be 
linear  (instead  of  quadratic).  The  approach  is  illustrated  in  figure  1. 

First,  some  concepts  and  notation  have  to  be  introduced.  Three  ingredients  are  funda¬ 
mental  to  the  construction:  The  B4zier  form  of  a  multivariate  polynomial,  the  identification 
of  the  appropriate  Bezier  nets,  and  the  identification  of  a  particular  parameter  with  each 
of  the  conditions  defining  the  interpolant. 

To  define  the  Bdzier  form  of  a  bivariate  polynomial  we  consider  a  general  triangle  with 
vertices  V,,  where  »  =  1,2,3.  A  general  point  P  in  is  expressed  in  terms  of  barycentric 
coordinates,  that  is 

3  3 

P=  where 

1=1  «=i 

Any  polynomial  q  of  degree  m,  say,  can  be  uniquely  expressed  in  its  Bizier  form  as 


t+j  +  fc=m 


(see  Farin,  1983).  We  will  have  to  differentiate  q.  Let  D  denote  any  (unnormalized) 
directional  derivative,  i.e.  D  =  -§1  where  e  is  a  vector  in  8#®.  Then 


Dq{P)  = 


»+y+fc=m— 1 


where 

Cijk  =  "»(ci+ij,fcZ?6i  +  c,,,+i,fcl>62  -f-  Ci,;,fc+iZ?63)  (2.2) 

(Subscripts  are  separated  by  commas  whenever  one  or  more  of  them  cannot  be  written  as 
just  one  character). 

The  Cijic  are  the  Bizier  ordinates  of  q.  They  can  be  associated  with  the  points 
Cijk  •'=  (>^i  +  jVz  +  kV3)/m,  The  points  (C,yfc,c,;*)  6  SR®are  said  to  be  the  control  points 
of  the  B4zier  polynomial,  and  the  set  of  control  points  constitutes  the  Bizier  net  of  q. 
If  two  polynomials  are  defined  on  two  neighboring  triangles,  then  the  resulting  piecewise 
function  is  continuous  across  the  interface  iff  the  Bezier  ordinates  on  the  interface  coincide. 
Thus  B4zier  ordinates  on  the  interface  can  be  identified  and  the  union  of  the  two  Bezier 
nets  can  be  considered  the  B4zier  net  of  the  piecewise  polynomial  function  defined  on  the 
union  of  the  two  triangles.  This  idea  will  be  expanded  below. 

We  now  proceed  in  6  steps: 

Step  0:  Forcing  Continuity  on  the  Macrotriangle.  (This  is  step  0  because  the  task  is 
accomplished  by  choosing  suitable  notation  rather  than  doing  analysis.) 

We  express  the  centroid  of  the  macrotriangle  by  V4  =  (Vi  +  1^2  +  V3)/3  and  denote 
(internal  and  external)  edges  by  Cij  :=  Vj  -  V,  (i.e.  c^y  is  the  line-segment  from  V,  to  Vj  ) 
where  i,j  =  1,2,3, 4.  The  microtriangle  where  6,  =  0  (i  =  1,2,3)  is  denoted  by  T,. 


Figure  1 

Stencil  for  Bivariate  Scheme 

•  position  given 
gradient  given 

LCBD  linear  cross  boundary 

derivative 


We  now  extend  the  notion  of  barycentric  coordinates.  Any  point  P  in  the  macrotri¬ 
angle  can  be  expressed  uniquely  as 

4  4  3 

P  =  where  II  *  “  1»2,3,4. 

1  =  1  i=l  «=1 

The  6,  can  be  thought  of  as  piecewise  linear  cardinal  functions  on  the  triangulation  of 
the  macrotriangle  into  microtriangles.  On  each  microtriangle,  at  least  one  of  6],  62,  or  63, 
will  be  zero,  and  the  other  coordinates  will  be  the  ordinary  barycentric  coordinates  from 
that  microtriangle  that  contains  P. 

We  will  have  to  differentiate  on  the  microtriangles  in  the  direction  of  (internal  and  ex¬ 
ternal)  edges.  Consider  a  general  triangle  spanned  by  the  set  { Va,  Vt,  Vc)  and  suppose  that 
{hjik}  C  {a,  6,  c}.  Then  it  follows  from  the  cardinal  properties  of  barycentric  coordinates 
that 

1^  =  Skj  -  Ski  (2.3) 

0€ij 

where  S  is  the  Kronecker  delta. 

The  piecewise  cubic  interpolant  (on  a  general  macrotriangle)  can  now  be  written  as 


«(f)  = 


E 

t+j+fe+/=3 


(2.4) 


where,  by  convention,  0°  =  1. 

The  resulting  generalized  B4zier  net  is  illustrated  in  figure  2.  The  function  q  defined 
in  (2.4)  is  a  continuous  piecewise  cubic  function  on  the  macrotriangle  with  19  parameters 

Note  that  internal  continuity  (i.e.  continuity  everywhere  on  the  macrotriangle)  has 
been  obtained  simply  by  the  judicious  choice  of  a  notation  that  makes  Bezier  control  points 
coincide  on  internal  edges. 

Step  1:  Forcing  Internal  Differentiability.  In  this  step,  we  impose  conditions  on  the 
Bezier  ordinates  that  will  ensure  that  the  piecewise  cubic  function  is  differentiable  ev¬ 
erywhere  on  the  macrotriangle.  The  standard  way  of  forcing  differentiability  between 
triangles  is  to  require  that  certain  groups  of  4  Bezier  control  points  each  be  colinear  (see 
Farin,  1983).  We  rederive  the  conditions  for  the  present  special  case  because  the  technique 
is  generally  useful,  and  because  it  is  needed  here  in  a  slightly  more  general  context  (i.e. 
differentiability  at  the  centroid). 

Consider  for  example  edge  €14.  Any  cross-boundary  derivative  on  that  edge  must  be 
independent  of  whether  it  be  evaluated  on  the  microtriangle  T2  or  the  microtriangle  T3.  It 
is  sufficient  to  consider  any  one  particular  cross-boundary  derivative,  because  any  other  one 
can  be  expressed  as  a  linear  combination  of  the  considered  one  and  the  tangential  derivative 
(in  the  direction  of  edge  614)  which  is  continuous  because  the  function  under  construction 
has  been  forced  to  be  continuous  everywhere  on  the  macrotriangle.  So  consider  for  example 
the  derivative  in  the  direction  of  edge  634. 


0210 


300 


Figure  2 

The  generalized  Bezier  Net 


In  order  to  obtain  an  appropriate  condition  we  carry  out  the  differentiation  on  both 
T2  and  Ts,  restrict  both  of  the  resulting  expressions  to  edge  en,  and  equate  the  coefficients 
of  the  resulting  quadratic. 

We  differentiate  on  directly  using  (2.1-3)  and  on  7*3  by  expressing  (34  =  -r  642, 
and  then  using  (2.1-3).  Carrying  out  this  plan  yields  the  equations 

C0003  =  (Ci002  C0IO2  CoOJ2)/3 

C1002  =  (ciioi  Cion  C2ooi)/3 
C200I  =  (C2100  +  C2010  C300o)/3 

Proceeding  similarly  on  the  other  two  internal  edges  yields  the  conditions: 

C0102  =  (cnoj  +  Com  co2oi)/3 
C020I  =  (ci200  +  C0300  +  C021o)/3 
C0012  =  (cion  +  Com  +  coo2i)/3 
C002I  =  (ci020  +  C0120  Co03o)/3 

Thus,  in  the  present  special  case,  the  differentiability  conditions  simply  require  that 
B4zier  ordinates  on  internal  edges  (in  the  interior  of  the  macrotriangle)  be  the  averages  of 
the  neighboring  Bezier  ordinates  in  the  directions  given  by  the  internal  edges!  Also  notice 
that  in  each  of  the  above  conditions  one  of  the  Bezier  ordinates  plays  a  special  role.  Each 
equation  can  be  thought  of  as  expressing  the  parameter  on  the  left  in  terms  of  the  right 
hand  side  which  involves  parameters  that  are  still  to  be  determined. 

At  this  stage,  q  is  a  piecewise  polynomial  function  that  has  12  parameters  still  available 
and  that  is  continuously  differentiable  on  the  macrotriangle. 

The  C’  conditions  have  a  remarkable  consequence:  At  the  centroid,  the  function  q  is 
actually  This  fact  was  first  observed  by  (Farin,  1983)  but  the  following  argument  will 
be  useful  in  the  trivariate  context:  Consider  for  example  the  microtriangles  T\  and  T3.  We 
have  to  show  that  three  independent  second  order  derivatives  are  continuous  when  going 
through  the  centroid  from  one  microtriangie  to  the  other.  So  consider  the  derivatives 


de24dei4'  de24^^34  de\4de34 


Of  these,  the  first  two  are  tangential  in  the  direction  of  634  which  is  shared  by  Ti  and 
T3.  They  are  hence  continuous.  Similarly,  the  third  derivative  is  continuous  between  the 
mierotriangles  T\  and  T2  since  it  is  tangential  in  the  direction  C34.  It  is  also  continuous 
between  T2  and  T3  since  it  is  tangential  in  the  direction  of  614.  (Obviously,  all  mixed 
partial  derivatives  commute).  Hence  it  is  continuous  between  Tj  and  T3.  This  completes 
the  argument. 

Step  2:  Interpolation  to  Vertex  Data.  The  key  observation  is  that  we  enforced  dif¬ 
ferentiability  at  the  vertices  of  the  macrotriangle  and  that  hence  we  need  to  interpolate 
to  each  datum  only  on  any  one  particular  microtriangie  that  is  convenient  for  that  da¬ 
tum.  With  this  observation,  interpolation  is  straightforward.  For  convenience,  we  think 
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of  the  data  as  being  given  in  terms  of  some  primitive  function  F.  Proceeding  similarly  as 
in  (Barnhill  and  Farin,  1981)  we  obtain: 


C3000  — 

F(P,) 

dF 

C2100  = 

Me, 2 

dF 

C2010  = 

^dexz 

Ctl300  = 

F(V2) 

dF 

Cl  200  = 

^de2i 

C02IO  = 

( 

5C23 

C0030  = 

F(V^3) 

dF 

C102O  = 

dF 

C0120  = 

dez2 

3C3noo)) 


/3 


3c3ooo))/3 


3co30o))/3 
+  3co3oo))/3 

+  3coo3o))/3 
+  3coo3o))/3 


At  this  stage,  we  have  obtained  an  interpolant  that  is  differentiable  on  the  macrotrian¬ 
gle.  The  construction  so  far  has  not  taken  into  account  smoothness  requirements  between 
macrotriangles,  and  there  are  still  three  parameters  at  our  disposal. 

Step  3:  Forcing  Global  Smoothness.  We  wish  to  apply  the  Clough-Tocher  scheme  on 
each  macrotriangle  in  a  given  triangulation  such  that  the  overall  function  is  differentiable 
everywhere  in  the  domain  covered  by  the  triangulation.  Along  an  external  edge  of  a 
macrotriangle,  the  interpolant  reduces  to  a  univariate  cubic  which  is  uniquely  determined 
by  the  four  data  at  the  endpoints  of  the  edge.  At  this  stage,  the  scheme  is  therefore  globally 
continuous  (but  not  differentiable). 

To  force  diflferentiability  it  is  again  sufficent  to  consider  one  particular  cross-boundary 
derivative  across  an  edge  joining  two  macrotriangles.  However,  the  only  meaningful  direc¬ 
tion  is  perpendicular  across  the  edge  because  that  edge  is  the  only  object  common  to  the 
triangles.  Since  the  interpolant  is  piecewise  cubic  the  cross-boundary  derivative  restricted 
to  the  edge  will  in  general  be  quadratic.  The  data  given  at  the  edge  impose  only  two 
conditions.  Therefore,  a  third  one  must  be  imposed  consistently  on  both  macrotriangles. 
It  is  natural  to  require  that  the  perpendicular  cross-boundary  derivative  be  linear.  This 
choice  maintains  the  maximum  degree  of  precision  (quadratic)  that  is  attainable. 

We  proceed  as  follows:  Let  [i,j,k)  be  a  cyclic  permutation  of  (1,2,3).  A  normal  to 
Cjk  can  then  be  written  as 

n,  =  ej4  -  itCjk  (2.5) 

where  7,  is  determined  by  the  requirement  that  n^ocjk  =  0.  i.e.  -y,  =  -  (c;40C;fc)/(cj*°Cjfc). 
We  then  differentiate  q  in  the  direction  of  n,  using  (2.5)  and  (2.1-3),  restrict  the  derivative 
to  edge  Cj*,  and  set  the  leading  coefficient  of  the  resulting  quadratic  equal  to  zero.  This 
yields  the  conditions 
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V 


Coni  =  ('>1  (-C0300  3co210  -  Scojjo  -r  Co03o) 

+  (-"C0300  2co210  ■“  C0120  -*■  C002I  +  C020i))/2 

Cion  =  ('72(“Co030  3ci020  ~  3C2010  C3000) 

+  ('C0030  +  2Ci020  ■"  C2010  C20OI  ■*"  C002i))/2 
Clioi  ==  ('y3(-C3000  3C2100  “  3ci200  Co30o) 

+  (-C3000  2C2100  “  C1200  -r  C200I  +  C020i))/2 

At  this  stage,  we  have  achieved  our  objective:  a  piecewise  cubic  polynomial  interpolant 
that  is  gobally  differentiable.  The  construction  used  up  all  parameters,  i.e.  the  interpolant 
is  unique.  In  the  trivariate  case,  however,  there  will  be  a  family  of  smooth  interpolants. 
In  that  case,  disposing  of  the  remaining  parameters  will  constitute  step  5. 

The  above  equations  for  the  coefficients  of  q  describe  a  linear  system.  However,  in 
this  case,  the  system  happens  to  be  triangular!  Thus  the  coefficients  can  be  computed  by 
Forward  Elimination,  e.g.  in  the  following  sequence:  C30001  C201O1  C2100,  co3oo»  C120O1  C0210, 
coo3o»  cio20)  C0120,  C2001,  co20if  coo2i»  Cnoi,  Cion,  Com,  C1002,  C0102,  C0012,  coooa- 

Each  of  the  parameters  is  naturally  associated  with  a  particular  condition,  and  we 
can  think  of  any  condition  as  eliminating  the  corresponding  parameter.  The  sequence  of 
evaluation,  and  the  type  of  condition  corresponding  to  a  parameter  is  illustrated  in  figure  3. 


3.  Construction  of  the  lYivariate  Scheme 

In  this  section,  we  follow  the  template  provided  by  the  preceding  one.  We  assume  that 
we  are  given  values  of  position,  gradient,  and  Hessians  at  scattered  points  in  SR®  (see  figure 
4.  for  the  stencil).  We  also  assume  that  the  domain  has  been  tesselated  into  tetrahedra, 
and  only  consider  the  interpolation  problem  on  a  general  tetrahedron. 

The  system  of  equations  determining  the  interpolant  has  been  derived  and  solved  using 
the  symbol  manipulation  language  REDUCE  (A.C.  Hearn,  1983).  Writing  down  all  the 
equations  and  their  solution  in  the  text  would  be  tedious  and  error  prone.  Instead,  the 
REDUCE  output  of  the  solution  has  been  reproduced  photographically  in  the  appendix. 
For  a  thorough  understanding  of  the  material  in  this  section,  the  reader  should  familiarize 
himself  with  the  organization  of  the  appendix,  and  keep  refering  to  it  when  reading  this 
section.  For  a  less  thorough  understanding,  this  section  may  be  read  by  itself. 


3.0  Notation  and  Internal  Continuity 

We  consider  a  general  macrotelrahedron  with  vertices  V^i,  ^2,  V3  and  V4,  and  denote 
the  centroid  by  =  (Fj  -+  V2  +  V'a  -t-  Vi][A.  The  location  of  a  general  point  P  is  expressed 

SIS 

5  54 

P  =  where  =  1,  ^>1  =  0,  and  6,  >0,  i  =  1,2, 3, 4, 5 

i=l  1=1  t=l 

and  the  function  to  be  constructed  is  of  the  form 


?(p)  = 


E 

«4j+fc+(+m  =  6 


5! 


m 

5 
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Figure  3 

Role  of  Coefficients  and 
Sequence  of  Evaluation 

•  interpolation 
o  internal  smoothness 
external  smoothness 


Note  that  our  interpolanl  will  be  piecewise  quintic.  Since  the  data  include  second 
order  derivatives,  on  any  edge  there  are  six  degrees  of  freedom  (position  and  first  and 
second  order  tangential  derivatives  at  each  vertex).  Thus,  using  any  degree  less  than  5 
would  imply  that  data  could  not  be  supplied  independently  at  the  vertices.  Also,  since  an 
arbitrary  number  of  edges  may  emanate  from  a  given  vertex,  the  specification  of  derivative 
data  —  which  usually  is  a  necessary  preprocessing  step  —  could  not  be  carried  out  on  an 
individual  tetrahedron,  or  a  set  of  points  close  to  the  given  vertex,  but  rather  has  to  take 
into  account  the  data  at  ail  other  points,  making  derivative  generation  a  process.  As  was 
pointed  out  in  the  introduction,  it  would  have  been  preferable  to  build  a  C*  scheme  for 
C*  data,  which  could  have  utilized  piecewise  cubic  functions,  but  attempts  to  do  so  have 
not  succeeded  to  date. 

At  this  stage,  q  is  continuous  on  the  macrotetrahedron,  and  contains  121  free  param¬ 
eters. 


3.1  Internal  Differentiability 

As  in  the  bivariate  case,  we  obtain  the  conditions  that  Bezier  ordinates  on  the  inter¬ 
faces  between  microtetrahedra  should  equal  the  average  of  their  three  neighbors  (in  the 
directions  of  the  internal  edges).  For  details  see  equations  no.  53-78,  91-111,  and  116-121 
in  the  appendix. 

At  this  stage,  q  is  differentiable  everywhere  on  the  macrotetrahedron,  and  contains 
121  -  53  =  68  free  parameters.  However,  there  is  a  subtle  difficulty:  We  contemplate  in¬ 
terpolation  to  vertex  data,  including  second  order  derivatives.  As  in  the  bivariate  case,  we 
would  like  to  interpolate  to  data  only  on  suitable  microtetrahedra,  and  enforce  interpola¬ 
tion  on  the  other  ones  by  smoothness.  Thus  we  have  to  enforce  second  order  differentiability 
at  the  vertices  of  the  macrotetrahedron,  although  the  scheme  under  construction  is  only 
CK 

Remarkably,  it  turns  out  that  the  scheme  is  already  C*  at  the  vertices,  in  spite  of  our 
not  having  made  any  effort  to  achieve  this.  The  crucial  fact  is  that  three  (rather  than  two) 
microtetrahedra  meet  at  each  vertex.  The  argument  is  similar  to  the  bivariate  case  where 
we  obtained  serendipitous  smoothness  at  the  centroid  of  the  macrotriangle: 

Consider,  for  example,  vertex  Vi  and  continuity  of  second  order  derivatives  between 
the  microtetrahedra  Jj  and  T3.  We  need  to  find  six  independent  second  order  derivatives 
that  are  continuous  at  Vj.  The  interface  between  T2  and  T3  is  the  triangle  spanned  by 
{k"i,V4,Vs}.  Because  the  scheme  is  already  C’,  all  second  order  derivatives  tangential  in 
that  triangle  will  be  continuous.  Thus  for  example 

a*  52  a2 

deu^ciA  ’  aei4aei5  ’  dei^dcis '  '  de^deis 

are  all  continuous  at  Vj.  We  need  one  more  derivative.  This  is  given  by 

a2 

de^^des^ 

which  is  continuous  because  it  is  tangential  between  T2  and  and  also  between  T3  and 
T4,  and  hence  between  T2  and  T4.  This  completes  the  argument. 


Tj,  ox  Figure  4 

iue  Stencil  for  the  t  • 

ior  the  Tnvariate  S( 
®  position,  gradient,  and 


It  is  noteworthy,  although  not  relevant  to  our  construction,  that  by  a  similar  procedure 
the  function  q  can  be  shown  to  be  at  the  centroid  (where  all  four  microtetrahedra  meet). 


3.2  Interpolation  to  Vertex  Data 

Having  forced  the  appropriate  degree  of  internal  smoothness,  interpolation  to  the 
vertex  data  is  again  straightforward.  40  parameters  are  required,  see  equations  1-40  in  the 
appendix.  At  this  stage,  the  scheme  does  interpolate  and  is  internally  C*.  We  have  to  use 
some  of  the  remaining  28  parameters  to  force  external  smoothness. 

3.3  Intertetrahedral  Smoothness 

External  smoothness  is  built  up  starting  at  the  edges  of  the  macrotetrahedron  and 
then  proceeding  to  its  faces. 

Continuity  across  edges.  This  has  already  been  achieved,  since  along  any  edge  the 
function  q  reduces  to  a  univariate  quintic  polynomial  which  is  determined  uniquely  by  the 
six  data  given  at  the  vertices  of  the  edge. 

Differentiability  across  edges.  Any  cross-boundary  derivative  across  an  edge  will  in 
general  be  univariate  quadratic.  However,  the  vertex  data  imply  only  four  conditions.  The 
resulting  ambiguity  is  removed  by  requiring  that  any  first  order  derivative  perpendicularly 
across  an  edge  be  in  fact  cubic.  This  condition  is  in  terms  of  the  edge  and  data  on  the 
edge  only,  and  nence  forces  differentiability  across  the  edges. 

We  exemplify  the  analysis  by  considering  edge  eij.  We  need  to  consider  derivatives  in 
two  independent  directions  perpendicularly  across  the  edge,  and  construct  one  direction 
in  each  of  the  two  microtetrahedra  joining  along  the  edge.  For  example,  consider  the 
microtetrahedron  T4.  We  express  a  normal  to  cu  as 


«123  =  ei3  +  7i23Cj2 


where 


Requiring  that  dq/dni23  be  cubic  yields  a  condition  for  022100-  Similarly,  considering 
ei2  a  part  of  T3  yields  a  condition  for  C22oio-  Ten  more  similar  conditions  are  obtained 
by  considering  the  other  five  edges  of  the  macrotetrahedron.  See  equations  41-52  in  the 
appendix  for  details. 

Continuity  across  faces.  This  is  implied  by  our  course  of  action  along  the  edges. 
On  each  face,  the  interpolant  reduces  to  a  specific  bivariate  quintic  interpolant  (Qig,  see 
Barnhill  and  Farin,  1981)  which  is  determined  uniquely  by  the  vertex  data. 

Differentiability  across  faces.  Consider  a  derivative  perpendicular  to  a  face.  In  general, 
this  will  be  a  bivariate  quartic  polynomial,  having  15  degrees  of  freedom.  Of  these,  9  have 
been  removed  by  vertex  interpolation  and  three  have  been  used  in  forcing  the  derivative 
to  be  cubic  along  edges.  The  last  three  degrees  of  freedom  are  removed  by  requiring  the 
derivative  on  the  face  to  equal  a  reduced  cubic  which  interpolates  to  the  nine  data  on  the 
face  (and  is  cubic  along  edges),  and  where  the  coefficient  corresponding  to  the  centroid  of 


»1123OCi2  =  0,  i.e.  7123  =  - 


gi3  o  C12 
C12  o  C12 
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f 


the  face  has  been  chosen  so  as  to  minimize  the  sum  of  the  squares  of  the  third  derivatives 
in  the  directions  of  the  edges. 

Specifically,  if  the  cross- boundary  derivative  on  a  face  spanned  by  V,,  Vj,  Vk  is  ex¬ 
pressed  as 

a  +  0-*-'i  =  3  ■  ■  '■ 


diw  =  (3(d210  +  <^201  "t"  di20  +  <^102  do21  +  <^012)  “  2(do03  +  ^fo30  +  d30o))/l2 

This  choice  maintains  quadratic  precision  of  the  derivative  (and  cubic  precision  of  the 
overall  interpolant  q). 

We  need  an  expression  for  a  direction  perpendicular  across  a  face.  Consider  for  ex¬ 
ample  face  1.  The  normal  is  expressed  as 

Til  =  625  +  0:1623  +  ^1634 

where  Oi  and  0i  are  determined  by  the  linear  equations 

Ml  0  623  =  fll  O  634  =  0. 

This  process  yields  three  conditions  on  each  face.  Detailed  expressions  are  listed  in 
equations  79-90  of  the  appendix. 

At  this  stage,  we  have  an  internally  and  externally  differentiable  interpolant  which, 
however,  still  contains  four  free  parameters. 

3.4  Condensation  of  Parameters 

We  still  have  to  dispose  of  the  four  parameters  coui2i  cioiiZi  Ciioi2»  and  611102-  In 
choosing  them,  it  is  desirable  to  meet  two  objectives; 

(1)  The  scheme  should  reproduce  exactly  polynomials  of  a  degree  as  high  as  possible  (i.e. 
cubics). 

(2)  The  interpolant  should  be  independent  of  a  relabeling  of  the  vertices. 

These  objectives  are  achieved  by  picking  the  remaining  parameters  so  as  to  minimize 
{Dfq{Vs))^  for  each  T,  where,  for  example 


D.  =  ( 


d  ^  ^  ^ 

de2i  ^635  ^645 , 


See  equations  112  through  115  in  the  appendix  for  the  detailed  conditions. 

3.5  Solution  of  the  Linear  System 

The  conditions  derived  in  the  preceding  subsections  constitute  a  linear  system  of  121 
equations.  In  the  bivariate  case,  the  corresponding  system  happened  to  be  triangular.  In 
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this  case,  it  turns  out  that  only  106  equations  can  be  solved  by  Forward  Substitution,  and 
that  the  remaining  15  equations  form  a  full  matrix.  The  overall  system  was  solved  using 
REDUCE,  and  the  results  are  listed  in  the  appendix. 


3.6  Precision  of  the  Trivariate  Scheme 

An  important  attribute  of  any  interpolation  scheme  is  its  degree  of  precision,  i.e.  the 
maximum  degree  of  a  polynomial  that  is  reproduced  exactly  by  the  scheme  if  the  data 
are  generated  by  the  polynomial.  Tracing  back  through  the  construction  it  is  easy  to  see 
that  our  scheme  is  precise  for  cubic  polynomials.  On  the  other  hand,  by  construction, 
derivatives  perpendicular  across  a  face  are  reduced  cubics  with  only  quadratic  precision. 
Therefore,  the  overall  scheme  cannot  have  more  than  cubic  precision.  A  parameter  count 
might  suggest  that  a  higher  degree  is  possible:  There  are  40  vertex  data,  and  a  trivariate 
quartic  has  only  35  coefficients.  However,  the  quartic  61626364  generates  homogeneous 
data  at  the  vertices.  Hence  any  linear  interpolation  process  will  yield  the  zero  function  as 
the  interpolant.  Thus  our  scheme  has  the  maximum  attainable  degree  of  precision. 


4.  Computational  Aspects 

In  this  section,  some  practical  aspects  are  discussed  as  they  have  been  incorporated 
into  a  (non-portable)  FORTRAN  research  code  that  can  be  obtained  from  the  author. 

1.  ConstructJou  of  the  Code.  The  present  scheme  shares  with  most  other  local  mul¬ 
tivariate  schemes  the  drawback  of  great  computational  complexity.  The  feasibility  of  pro¬ 
gramming  such  schemes  from  printed  formulas  is  called  into  question.  Instrumental  in  the 
writing  of  the  code  was  a  feature  of  the  symbol  manipulation  language  REDUCE  that 
allows  it  to  write  formulas  in  FORTRAN  syntax  after  deriving  them  from  much  simpler 
ingredients. 

2.  Labeling  of  coefficients.  Having  coefficients  with  five  subscripts  is  convenient  for 
the  analysis  but  not  suitable  for  a  computer  Implementation.  One  would  like  to  use  a 
linearly  indexed  array  for  storing  the  coefficients.  We  address  the  problem  somewhat  more 
generally: 

Let  I  denote  the  set  of  all  Jb-tuples  of  non-negative  integers  that  add  up  to  m.  There 
are  M  :=  *1*)  elements  in  7.  A  bijective  correspondence  between  7  and  {1,2,  ••  • ,  M} 

is  defined  in  terms  of  binomial  coefficients  by: 


N(i 


i>  ■ 


ik)  =  1  + 


^U-i  ^  -f  +  *fc-i  +  ^  ^  ^*1 


fc  —  1 


) 


For  the  tetrahedral  Clough-Tocher  scheme,  m  =  A:  =  5,  Af  =  126,  and  the  indices 
iV(l,l,l,l,2)  =  50,  iV{l,l, 1,2,1)  =  99,  Ar(l, 1,2,1, 1)  =  98,  Ar(l, 2, 1, 1, 1)  =  95,  and 
iV(2, 1, 1, 1, 1)  =  85  are  unused. 

3.  Identification  of  the  microtetrahedron.  Suppose  the  macrotetrahedron  in  which 
the  point  of  evaluation,  P,  say,  resides  has  been  identified.  (This  is  done  by  evaluating 
the  barycentric  coordinates  with  respect  to  the  tetrahedron  and  checking  that  they  are 
all  non-negative.)  Now  the  microtetrahedron  containing  P  must  be  identified.  Denote 
the  barycentric  coordinates  by  6,,  bj,  6*,  6j  where  {i,j,kj}  -  {1,2, 3, 4}  and  assume 


bi  =  min{bi,bj,bk,bi}.  Then  P  is  contained  in  T,.  Denoting  the  barycentric  coordinates 
of  P  with  respect  to  T,  by  ds,  dj,  dk^  di  it  is  easy  to  see  that  df,  =  46,,  dj  =  bj  -  6,, 
dk  =  bk  -  6,,  di  —  bi  -  bi. 

4.  Recursive  Evluation  of  Interpolant.  The  Bezier  form  of  a  polynomial  allows  for 
a  particularly  simple  recursive  form  of  the  evaluation.  For  details  see  Farin,  1983. 

5.  Efficient  evaiuation  of  coefficients.  Coefficients  are  evaluated  only  on  macrotetra- 
hedra  where  they  are  needed.  The  user  has  the  option  of  storing  them  after  evaluation,  or 
reevaluating  them  when  needed. 

6.  Verification  of  Code.  The  pilot  code  developed  by  the  author  has  been  tested  by 
the  analyzing  tool  MICROSCOPE  (Alfeld  and  Harris,  1984)  and  it  was  verified  that  the 
code  does  indeed  possess  the  smoothness,  interpolation,  and  precision  properties  that  are 
implied  by  the  mathematical  construction. 

Conclusions 

The  scheme  developed  here  is  the  first  explicitly  given  piecewise  polynomial  C*  inter- 
polant  for  tetrahedral  data.  It  has  the  drawback  of  requiring  data. 
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Appendix:  Coefficients  of  Interpolant 

In  this  appendix  we  list  the  formulas  as  they  have  been  computed  by  REDUCE. 
Partial  derivatives  of  the  Function  F  have  been  denoted  by  subscripts,  for  example 

F.4(V,)= 

oei4 


d^F 

0€i^0t\4 


1:  c 50000  = 

2:  C40010  = 
3:  C40100  = 
4:  C41000  = 
5-'  C30020  = 
6:  C30200  = 
7;  C32000  = 
8:  C30UO  = 
9:  C31010  = 
10:  C31100  = 
11:  cosooo  = 
12:  C04010  = 
13:  C04100  = 
14:  c 14000  = 
15:  C03020  = 
16:  C03200  = 

17:  C23000  = 
18:  C03110  = 
19:  C13010  = 

20:  C13100  = 
21:  C00500  = 
22:  C00410  = 
23:  C0I400  = 
24:  C10400  = 
25:  C00320  = 
26:  C02300  = 


(5C5oOOO  +  ^’i4{V'i))/5 

{5C50000  +  P’i3(Vi))/5 

(5C50000  +  Fi2(V'i))/5 

(-2OC500OO  +  4OC400IO  +  f^I4I4{l^l))/20 

(-2OC5000O  +  4OC40100  +  f^l313(l^l))/20 

(-2OC50000  +  4OC41000  +  ^’i212(1^i))/20 

(-2OC500OO  +  2OC4010O  +  2OC400IO  +  ■^^1413(1^i))/20 

(-2OC50000  +  2OC41000  +  2OC40010  +  ^mi2(1^i))/20 

(-2OC500OO  +  2OC4100O  +  20C40100  +  i^l312(l^l))/20 

F(V2) 

(5Co5000  +  f^24(1^2))/5 

(5Co5000  +  F23(V’2))/5 

(5C0S000  +  F2i(V2))/5 

(-20co5OOO  +  40Co4010  +  F2424(1^2))/20 

{-20Co6000  +  40Co4100  +  F2323(1^2))/20 

(40Ci4000  -  20Co5000  +  F212i(F2))/20 

(-20Co6000  +  20co4IOO  +  20co4010  +  F2423(1^2))/20 

(20Ci4000  -  20co5000  +  20Co4010  +  F242|(F2))/20 

(20c  14000  -  20co5000  +  20Co4100  +  F232i(V'2))/20 

F(V3) 

(5Coo500  +  F34(V3))/5 

(5Coo500  +  F32(V3))/5 

(5Co0800  +  F3i(V3))/5 

(-20COOSOO  +  40Co04IO  +  F3434(V3))/20 

(40Co1400  ~  20Co0500  +  F3232(1^3))/20 


27:  C20300  =  (40Cjo4OO  -  20CoO5OO  +  i^313l(^3))/20 

28:  C01310  =  (20couoo  —  20coo50o  +  20coo4io  ^3432(V'3))/'20 

29:  C10310  =  (20cio4oo  -  20coo50o  +  20coo4io  +  ^343i{V^3))/20 

30:  C11300  =  (20Ci04OCi  +  20coi400  -  20C005OO  ^323l(V'3))/20 

31:  coooso  =  F{V4) 

32:  C00140  =  (Scoooso  +  ^43(^^4))/5 

33:  C01040  =  (Scoooso  +  ^42('^4))/5 

34;  c  10040  =  (Scoooso  +  ^^41  (^4)) /5 

35:  C00230  =  (40Coo14O  ~  20coo050  +  i^4343(^4))/20 

36:  C02030  =  (40Coi040  ~  20coo050  +  ^4242(^4))/20 

37:  C20030  =  (40cioo40  -  20coooso  +  ^4i4i(^4))/20 

38:  C01130  =  (20CoiO4O  +  20cooi40  -  20coooso  +  ^4342{^4))/20 

39:  C10130  =  (20cioo40  +  20cooi4o  -  20coooso  -*•  ■f434i(V4))/20 

40:  C11030  =  (20cioo40  +  20coio4O  ~  20coooso  +  ^424i(^4))/20 

41:  C22IOO  =  ('yi23(C50000  -  5C41000  +  IOC32OOO  “  10C23000  +  5Ci4000  -  Cosooo)  +  C50000  - 

4C41000  -  «40100  +  6C32000  +  4C31100  -  4C23000  +  C14000  +  4Ci3)00  “  C0410o)/6 

42:  C23010  =  ('1124(^50000  -  5C41000  +  IOC320OO  —  IOC23OOO  +  Scuooo  ~  cosooo)  +  C50000  - 
4C41000  -  C40010  +  6C32000  +  4C31010  -  4C23000  +  C14000  +  4Ci30lO  -  C0401o)/6 

43:  C21200  =  ('yi32(C50000  -  5C40100  +  IOC302OO  —  10C20300  +  5Ci0400  ~  C00600)  +  C50000  - 
C4IOOO  ~  4C40100  +  4C3iioo  +  6C30200  ~  4C20300  +  4Cii300  +  C10400  “  Coi40o)/6 

44;  C20210  =  (7134(c50000  -  5C40100  +  IOC302OO  -  IOC2030O  +  5ci0400  -  Coosoo)  +  C50000  - 
4C40100  —  C40010  +  6C30200  +  4C30110  ~  4C20300  +  C 10400  +  4c  10310  -  Co041o)/6 

45:  C21020  =  ('7142(C50000  -  5C40010  +  IOC3002O  -  IOC20030  +  5ci0040  ~  CoOOSo)  +  C50000  - 
C4IOOO  “  4C40010  +  4C310I0  +  6C30020  -  4C20030  +  4Cii030  +  C10040  -  Coi04o)/6 
46:  C20130  =  ('7i43(c60ooo  -  5C40010  +  IOC30020  -  10c2oo30  +  5c  10040  -  Coooso)  +  C50000  - 
C40100  -  4C40010  +  4C30110  +  6C30020  -  4C20030  +  4Cioi30  +  C10040  -  Coo14o)/6 
47:  C12200  =  ('1231  (cosooo  -  5C04100  +  IOC032OO  -  IOC0230O  +  5Coi400  -  CooSOo)  -  Ci4000  + 
4Ci3ioo  +  4Cii300  -  C10400  +  Cosooo  -  4co4100  +  6co3200  -  4Co2300  +  Co140o)/6 

48:  C02210  =  ('y234(co500o  -  5co410o  +  10co320o  -  10co23oo  +  ScoMOO  -  Coosoo)  +  cosooo  - 
4Co4100  -  C04010  +  6Co3200  +  4co3110  -  4Co2300  +  Couoo  4-  4Coi310  —  Co041o)/8 

49:  C12020  =  ('7241  (cosooo  -  5co40io  +  IOC03020  -  IOC02030  +  5coio40  -  coooso)  -  Cuooo  + 
4Ci3010  +  4Cii030  “  C10040  +  Cosooo  -  4Co4010  +  6co3020  ~  4Co2030  +  Coi04o)/6 

50:  C02120  =  ('7243(co5ooo  -  5co40io  +  10co3020  -  10co2030  +  5coio40  -  Coooso)  +  Cosooo  - 
C04100  -  4Co40IO  +  4Co3IIO  +  6co3020  -  4Co2030  +  4Coi130  +  C01040  -  Coo14o)/6 
51:  C10220  =  ('734i(C00S00  -  5co0410  +  10coo320  -  10CoO23O  +  5coo140  -  Coooso)  ~  C10400  + 
4Ci0310  +  4Cio130  -  CI0040  +  Coosoo  -  4co0410  +  6CoO320  -  4coo230  +  Coo14o)/6 
52:  C01220  =  (7342(coosoo  -  5coo4io  +  10coo32o  -  10coo23o  +  5cooi40  -  Coooso)  -  cqhoo  + 
4coi31o  +  4coii30  -  cok)4o  +  Coosoo  -  4coo4io  +  6C00320  ~  4coo230  +  cooi4o)/6 


53;  C4000I  =  (csoooo  +  C41000  +  C40100  +  C4ooio)/4 
54:  C31001  =  (c4iooo  +  C32000  +  caiioo  +  C3ioio)/4 
55:  C22OOI  ==  (C32OOO  +  C23000  "t-  C22IOO  +  C2201o)/4 
56:  Ci30Cil  =  (C23000  +  C14000  +  ClSlOO  +  C1301o)/4 
57:  C04001  =  (ci4ooo  +  cosooo  +  C04100  +  co40io)/4 
58:  C30U11  —  (C40100  +  C3noo  +  C30200  ■+  C3oiio)/4 
59:  C03IOI  =  (Cl3100  +  C04IOO  +  C03200  +  co3iio)/4 
60:  C20201  =  (C30200  +  C21200  +  C20300  +  C2021o)/4 
61:  C022OI  =  {ci2200  +  C03200  +  C02300  +  Co221o)/4 
62:  C10301  =  (c20300  +  cnsoo  +  C10400  +  cio3io)/4 
63:  C01301  =  (cn300  +  C02300  +  C01400  +  coi3io)/4 
64:  C004OI  =  (c 10400  +  C01400  +  cqosoo  +  coo4io)/4 
65:  C30011  =  (C40010  +  C31010  +  C30110  +  C3002o)/4 
66;  C03011  =  (ci30io  -t  C04010  +  C03110  +  co3o2o)/4 
67:  C003U  =  (cio3io  +  coiaio  +  coo4io  +  coo32o)/4 
68:  C20021  =  (C30020  +  C21020  +  C20120  +  C2003o)/4 
69;  C02021  =  (c  12020  +  C03020  +  C02120  +  C0203o)/4 

70:  C00221  =  (ci0220  +  C01220  +  C00320  +  Coo23o)/4 

71:  C10031  =  (C20030  +  C11030  +  C10130  +  Ci004o)/4 
72:  C01031  =  (ciio3o  +  C02030  +  conao  +  coio4o)/4 
73:  C00131  =  (cioi3o  +  Couao  +  C00230  +  cooi4o)/4 

74:  C0004I  =  (Cl0040  +  C01040  +  C00140  +  cooo5o)/4 

75:  C30002  =  (^40001  +  C31001  +  C30101  +  C3ooii)/4 

76:  C03002  =  (Cl3001  +  C04OOI  +  C031OI  +  Co301l)/4 

77:  C00302  =  (c  10301  +  C01301  +  C00401  +  coo3ii)/4 

78:  C00032  =  (cioo3i  +  Coioai  +  cooiai  +  Cooo4i)/4 

79:  C02111  =  (oi(23co5ooo~61co41oo~38co4010  +  59co320o+86co3iio  +  21co3020~43co2300~ 
48co2210  -  2IC02120  ~  22Co2030  +  31Co1400  +  2Co1310  -  15Coi220  +  24Coi130  +  9Coi040  - 
9Co0500  —  2Co0410  +  15co0320  “  2Co0230  —  9Co014o)  +  0i  (23co4100  ~  23Co4010  —  38Co3200  + 
38co3020  +  21Co2300  +  27co2210  ~  27co2120  ~  21Co2030  ~  22coi400  22co1310  —  22Coi130  + 

22co1040  +  9co0500  -  7Co0410  -  17co0320  +  17Co0230  +  7co0140  -  9co005o)  +  23co5000  - 
38Co4100  ~  38Co4010  ~  23C04001  +  21C03200  +  48Co31  10  +  38Co310I  21Co3020  +  38Co3011  “ 

22co2300  —  2IC022OI  ~  22Co2030  "  21Co2021  9co1400  +  2Co1310  22co1301  ~  15Co1220  + 

2Coi130  +  9coi040  +  22Coi031  -  9Co0401  -  2Co0311  +  15Co0221  ~  2coo131  “  9C0004i)/48 
80:  C01211  =  (oi(9co5ooo  “  31co4ioo  2co4oio  +  43co320fi  ~  2co3i  10  ~  15co3020  “  59co230o  + 
48Co2210  +  15Co2I20  +  2Co2030  +  61Coi400  ■"  86Coi310  21Coi220  -  24coii30  9Coi040  — 

23Co0500  +  38Co0410  -  21Co0320  +  22C0023O  "  9coo14o)  +  ^1  (9co4100  ~  9C04010  ~  22Co3200  + 
24Co3I10  ~  2co3020  +  21Co2300  -  21Co2210  -  15Co2120  15C0203n  -  38Co1400  +  86Coi310  ~ 


-20- 


48coi220  2con3o  ~  2coi04Ci  23cooioo  ~  6lcoo4Jrt  SQcnoszo  ~  43coo230  +  31cooi40  ~ 

9Co006o)  9c05000  ”  22C04Ki0  +  2co40JO  ~  9c04001  21Co3200  22co3101  —  15Co3020  ~ 

2co30U  “  38co230Cj  +  48co2210  ~  21Co2201  +  2co2030  15co2021  23Coi400  “  38Co1310  + 

38coi3oi  21coi22o  ~  22coii3o  +  9coio40  “  2coio3i  ~  23coo40j  38coo31i  “  21coo22j 
22cooi3i  ~  9tooo4i)/48 

81:  CCI112J  =  (oi  (9to5ooi;>  ~  7co4iof' ~  22cci40io ■■  17co3200  22co31i(i 21co3020  17co230o  + 

27co212(.'  ~  38co2ci3o  +  7coi40o  ~  22eoi3io  ~  27c(ii22ci  23coio40  ~  9coo50o  +  22coo4io  ~ 

21Coo320  38Co0230  “  23co014o)  +  /3j(9co4100  “  9Co4010  2Co3200  ~  24Co3110  22C03020  ~ 

15Co2300  +  15C02210  +  21Co2120  -  21co2030  +  2co1400  ~  2Coi310  +  48^01220  ~  86coii30 
38Coi04a  •i'  9Co0500  "  31C00410  +  43co0320  ~  59Co0230  +  61^00140  ~  23CoOOSo)  9co5000 
2c04100  ~  22C04010  ~  9co4001  ~  15Co3200  —  2Co3101  +  21C03020  22C03011  +  2C02300  + 

15Co2201  +  48Co2120  ~  38co2030  ~  21Co2021  +  9co1400  ~  22coi310  ~  2Co1301  +  21coi220  ~ 
38coii3o  +  23coio40  +  38coio3i  "9600401  +  22000311  ~21coo22i  +38cooi31  ~ 23cooo4i)/48 

82:  C20111  =  (O2(23C40100~23C40010“38C30200+38C30020  +  21C20300  +  27c20210"27C20120~ 
21C20030  “  22Ci0400  +  22Cio310  “  22Cio130  +  22Cio040  +  9Co0500  ~  7Co0410  ~  17Coo320  + 
17co0230  +  7coo140-9coo05o)+/J2(~23C60000  +  38C40100  +  61C40010~21C30200~86C3oiio~ 
59C30020  +  22620300  +  21620210  +  48620120  +  43620030  ~  9eio400  ~  246x0310  +  156x0220  ~ 
2cio130  ~  316X0040  +  9600410  +  2600320  ~  15600230  +  2eo0140  +  9Co006o)  +  23640100  ~ 
23640001  "  38630200  ~  38630110  +  38630X01  +  38630011  +  21620300  +  48620210  “  21620201  + 
21620120  ~  21620021  ~  226x0400  +  226x0301  ~  226x0130  +  226x0031  +  QoqOSOO  +  2600410  “ 
9600401  ~  15600320  ~  2600311  +  2600230  +  15600221  +  9600140  “  2600131  ~  9eoo04l)/48 

83:  6X0211  =  (0!2(9C40100  ~  9640010  -  22630200  +  24630110  ~  2630020  +  21620300  “  21620210  “ 
15620120  +  15620030  “  386x0400  +  866x0310  “  486x0220  +  26x0130  ~  26x0040  +  23600500  ~ 
61600410+  59600320"  43600230 +  31600140 “96 O005o)  +  /?2(~96 50000  +  22640100+ 7640010“ 
21630200  "  226301X0+  17630020  +  38620300  ~  27620210  "  17620030  "  236x0400  +  276x0220  + 
226  10130  "  76X0040  +  23600410  "  38600320  +  2l6oo230  "  22600140  +  96o005o)  +  9640X00  " 
9640001  "  22630200  +  2630110  +  22630101  "  2630011  +  21620300  "  21620201  "  15620120  + 
15620021  "  386x0400  +  486x0310  +  386x0301  +  26X0130  "  26x0031  +  23600500  "  38600410  " 
23600401  +  21600320  +  38600311  "  22600230  ~  2l6oo221  +  9600140  +  22600131  ~  9eo004l)/48 

84:  610121  =  («2(964oioo  "  9640010  +  2630200  "  246301 10  +  22630020  "  15620300  +  15620210  + 
21620120  "  21620030  +  26i0400  "  26i0310  +  48Ci0220  "  80610130  +  38C|0040  +  9600500  " 
31600410  +  43co0320  "  59600230  +  61600140  "  236oo05o)  +  ^2( "9650000  "  2C40100  + 
31640010  +  15630200  +  2630110  "  43630020  ~  2620300  "  15620210  "  48620120  +  59620030  " 
9610400  +  24610310  "  21610220  +  80610130  ~  61Ci0040  +  9Co04IO  "  22600320  +  2l6o0230  " 
38600140  +  23co005o)  +  9640100  "  9640001  +  2630200  "  22630110  "  2630101  +  22630011  " 
15620300  +  15620201  +  21620120  "  21620021  +  2eio400  ~  2610301  +  486i0220  "  386ioi3o  + 
38610031  +  9600500  "  22600410  -  9600401  +  21co0320  22co0311  "  38600230  "  21600221  + 

23600140  +  38600131  "  2360004i)/48 

85:  621011  =  (o3("23650000  +  38641000  +  6IC40010  "  2IC32000  "  80631010  "  59630020  + 
22623000  +  2I622010  +  48621020  +  43620030  "  96i4000  ~  246x3010  +  15612020  ~  2Ciio30  " 
31610040  +  9604010  +  2603020  "  15co2030  +  2601040  +  9eoo05o)  +  /?3(23e50000  "  6IC41000  — 
38640010  +  59632000  80631010  +  21630020  ~  43623000  ~  48622010  ”  21621020  "  22620030  + 

31  ^14000  +  26x3010  "  156x2020  +  246ii030  +  9<'l0040  ~  9co5000  "  2604010  +  15Co3020  " 


2co203o  -  9coio4o)  +  23c4ocao  -  23c4oooi  -  38031010  ^  38c3iooi  -  38c30020  +  38c3ooii 
2IC22OIO  “  2IC2200I  +  48C2  102(1  2IC2003O  21020021  “  220i3010  +  22c  13001  "•  22c  10040  + 

22c 10031  T  9C04010  “  9co4001  ^  2co3020  “  2co3011  -  15c02030  +  13co2021  +  2CoI040  “ 
2Coi031  +  9Coo050  “  9Co004i)/48 

86:  C1201I  =  (Q;3( -Qcsoooo"^  22C4100i,i*  "C40(pK'- 2IC32000  -  22031010+ 17C30020  + 38C23000~ 
27022010  “  17C20030  “  23Ci40oo  27ci2ii2(i  +  22ciici3('  —  7oio040  +  23co4010  “  38Co3020 
21Co2030  -  22Coi040  +  9co005o)  y33(9C50000  -  31041000  2c40010  +  43032000  —  2031010  ~ 

15^30020  ~  39023000  +  48022010  +  15021020  +  2020030  +  61^14000  ~  86^13010  +  21Ci2020  ~ 
24Cii030  +  9Cio040  ~  23Co5000  +  38Co4010  “  21C03020  t-  22co2030  ~  9Coi04o)  +  9C40010  ~ 
9C40001  “  22031010  +  22031001  +  2030020  “  2030011  2IC22010  “  2IC22OOI  —  15C20030  + 

15020021  —  38ci3010  +  38ci3001  +  48ci2020  +  2ci0040  ~  2Ci0031  +  23Co4010  ~  23Co4001  ~ 
38Co3020  +  38Co3011  +  21Co2030  ■“  21Co2021  ~  22Coi040  22co1031  +  9Coo050  ~  9C0004i)/48 

87:  C1102I  ==  (oi3(  — 9050000  ~  2041000  +  3IC40010  +  15C32000  +  2031010  “  43030020  ~  2023000  ~ 
15022010  ~  48021020  +  59020030  ~  9Ci4000  +  24Ci3010  “  21Ci2020  +  86Cii030  “  61Ci0040  + 
9co4010  -  22Co3020  +  21Co2030  -  38coi040  +  23coo05o)  +  /S3(9c 50000  ~  7C4iooO  “  22040010  ~ 
17C32000  +  22031010  +  2IC30020  +  17023000  +  27021020  ”  38C20030  7ci4000  ~  22Ci3010  ~ 

27Ci2020  +  23ci0040  ~  9coSOOO  +  22Co4010  ~  21C03020  +  38co2030  ~  23coi04o)  +  9040010  ~ 
9C40001  +  2031010  —  2031001  ~  22C30020  +  22C30011  “  15022010  +  15022001  +  21020030  “ 
2I020021  +  2ci3oio  -  2ci3ooi  +  48ciio30  -  38cioo40  +  38cioo3i  +  9co40io  —  9co400i  " 
22co3020  +  22co3oii  +  21 002030  “21co2021  —  38coio40  +  38coio31  +  23000050“  23cooo4i)/ 48 

88:  C21101  =  (o!4(23c5oooo“61<^4iooo“38c4oioo+59c32ooo+86c3iioo  +  21c30200“43c230oo“ 
48022100  “  2IC21200  “  22020300  +  31Ci4000  +  2Ci3100  “  15Ci2200  +  24Cii300  +  9Ci0400  “ 
9co8000  “  2Co4100  +  15Co3200  “  2Co2300  “  9co140o)  +  /34(23c41000  “  23040100  “  38032000  + 
38030200  +  21023000  +  27022100  ”  27021200  “  21020300  “  22Ci4000  +  22Ci3100  “  22cii300  + 
22ci0400  +  9co5000  “  7Co4100  “  17C03200  +  17co2300  +  7coi400  “  9Co050o)  +  23C50000  “ 
38041000  “  38c4oioo  “  23c4oooi  +  2IC32000  +  48c3i  100  38c3iooi  +  2IC30200  +  38c3oioi  ~ 
22023000  “  2IC2200I  “  22020300  “  2IC2020I  +  9ci4000  +  2Ci3I00  +  22Ci300I  “  15Ci2200  + 
2Cii300  +  9Ci0400  +  22Ci0301  “  9C0400I  “  2Co3101  +  15Co2201  “  2Co1301  “  9C0040i)/48 

89:  012101  =  (0:4(9050000  “  3IC41000  +  2C40100  +  43C32000  “  2031100  “  15C30200  “  59C23000  + 
48C22100  +  15C21200  +  2C20300  +  61Ci4000  “  86ci3100  +  21Ci2200  “  24Cii300  +  9Ci0400  ~ 
23Co5000  +  38co4ioo  “  21C03200  +  22Co2300  “  9Co140o)  +  /?4  (9041000  “  9C40100  “  22032000  + 
24031100  “  2030200  +  21023000  “  2IC22100  “  15021200  +  15020300  “  38Ci4000  +  86Ci3100  “ 
48Ci2200  +  2Cii300  —  2Ci0400  +  23C05000  “  61co4100  +  59Co3200  “  43Co2300  +  31Co1400  ~ 
9Co050o)  +  9050000  “  22C41000  +  2040100  “  9C40001  +  2IC32000  +  22C31001  “  15C30200  “ 
2030101  “  38023000  +  48022100  “  2IC22OOI  +  2020300  +  15020201  +  23Ci4000  “  38Ci3100  + 
38ci3001  +  21ci2200  “  22Cii300  +  9Ci0400  ~  2Cio301  “  23Co4001  +  38Co3101  “  21Co2201  + 
22coi30i  “  9coo40i)/48 

90:  C112OI  =  (0:4(9050000  “  7041000  “  22C40100  ~  17032000  +  22C3i  100  +  21030200  +  17023000  + 
27C21200  “  38020300  +  7ci4000  “  22ci3100  ~  27oi2200  +  23ci0400  “  9Co5000  +  22co4100  “ 
21Co3200+  38Co2300  “  23Co140o)  +  ^^4  (9041000  “  9040100  +  2032000  “  24031100  +  22C30200  ~ 
15023000  +  15022100  +  21021200  “  21020300  +  2Ci4000  “  2Ci3ioo  +  48Ci2200  “  86Cii300  + 
38Ci0400  +  9C0S000  “  31Co4100  +  43Co3200  “  59Co2300  +  61Co1400  “  23Co050o)  4  Qcsoooo  + 
2041  000  “  22040100  “  9040001  -  15032000  “  2031001  +  2IC30200  +  22030101  +  2023000  + 


15C22001  +  48C21200  ~  38C20300  ~  2IC202OI  9Ci4000  ■"  22Ci3100  “  2Ci3001  +  21Ci2200  ~ 
38c  ij 300  23c  10400"^  38c  10301  —  9c 04001  22co3ioi  ~  21c 02201  +  38coi30i  ~  23coo40i)/48 

91:  C21002  =  (C31001  +  C22001  +  C2noi  -r  C2ioii)/4 
92:  C12OO2  =  (c22001  1-  C13001  C12IOI  Ci201l)/4 

93:  C20102  =  (C30101  ->■  C2noi  C20201  C2oiii)/4 

94:  002102  =  (Cl2101  Cosioi  C022OI  C0211l)/4 

95:  C10202  =  (C20201  +  cn2oi  +  C10301  cio2n)/4 

96:  C01202  =  (cii20i  +  C02201  -i-  C01301  -r  coi2ii)/4 

97:  C20OI2  =  (C3001I  +  C21011  +  C2OIII  -t-  C2002l)/4 

98:  C02OI2  =  (Cl2011  +  C030U  +  C02IU  C0202l)/4 

99:  C00212  =  (ci02ii  +  C01211  +  coo3n  t  coo22i)/4 

100:  C10022  =  (c2002i  +  Cl  1021  +  C10121  +  cioo3i)/4 

101:  C01022  =  (cuo21  +  C02021  +  C0112I  +  C0103i)/4 

102:  C00122  =  (ci0121  +  C0112I  +  C0022I  +  C0O13i)/4 

103:  C20003  =  (C30002  +  C21002  +  C20)02  +  C20012)/4 

104:  C02003  =  (ci2002  +  C03002  +  C02102  +  C02012)/4 

105:  C00203  =  (ci0202  +  C01202  +  C00302  +  Co0212)/4 

106:  C00023  =  (ci0022  +  C01022  +  C00122  +  Coo032)/4 

107:  C|ooo4  =  (255c4oooi  +  68OC31001  +  68OC30101  +  68OC30011  “  3O6OC30002  +  IO2OC22001  + 
IO2OC21101  +  1020c2iou  ”  4032c2ioo2  +  1020c2O2oi  +  1020c2oiii  “  4032c2oio2  + 
IO2OC20021  “  4032c2ooi2  +  8838c2ooo3  +  680ci3ooi  +  1020ci2ioi  +  1020ci2oii  ~ 
4032ci2oo2  +  1020ch2oi  +  1020ciio2i  +  680cio3oi  +  1020cio2ii  ~  4032cio202  + 

1020cio]2i  +680cioo31  ~4032cioo22  +  251co40oi  +664co3]oi  +664co30ii  —  3012co3002  + 
996co220i  +  972co2ih  ~  4032co2io2  +  996co202i  ~  4032co20i2  +  8478co2003  +  664coi30i  + 
972coi2ii  ~  4032coi202  +972coii2i  +  664coio3i  ~  4032coio22  +  251coo40i  +  664coo3ii  ” 
3012coo3O2  +  996coo221  ~  4032coo2i2  +  8478coo203  +  664cooi3i  ~  4032cooi22  +  251cooo4i  “ 
3012cooo32  +  8478cooo23)/1008 

108:  C01004  ==  (— 392c4oooi  “  68OC31001  ~  I228C30101  ~  I228C30011 +  4704c3ooo2“  IO2OC22001  ” 
IO2OC21101  ~  1020c2ioii  +  4032c2ioo2  ~  1842c2020i  ~  2664c2oui  +  4032c2oi02  “ 
1842c2oo2i  +  4032c2ooi2  ~  21168c2ooo3  —  68OC13001  ~  1020ci2ioi  —  1020ci2oii  + 
4032ci2oo2  ■"  1020cii2oi  ~  1020ciio2i  ~  1228cio30i  ~  2664cjo2ii  +  4032cio202  ~ 
2664cioi2i  ~  1228cioo3i  +  4032cioo22  "i"  35532ciooo4  ~  114co400i  ~  II6C03101  ~ 
116co30i  1  +  1368co3oo2  ~  174co220i  +  672co2ii  1  +  4032co2io2  ~  174co202i  +  4032co2oi2  + 
3852co2oo3  ~  116coi3oi  +  672coi2n  +  4032coi202  +  672coi  121  ~  1 16coio3i  +  4032coio22  ~ 
251coo401  ~  664coo3U  +  3012coo3n2  ~  996coo22i  +  4032coo2i2  ~  8478coo203  ~  664cooi3i  + 
4032cooi22  ~  251cooo41  +  3012cooo32  ~  8478cooo23)/34524 
109:  C00104  =  (~392c4oooi  ~  1792c3iooi“  116c3oioi~  1228c3ooii +  4704c3ooo2“ 2688022001  “ 
1020c2noi  —4356c  2 101 1  +  4032c2ioo2  —  174c2020i  +672c2oii  1  +  4032c  20102  —  1842c2oo2i  + 
4032c2ooi2  —  2II68C20003  —  1792ci3ooi  —  1020ci2ioi  —  4356ci2oii  +  4032ci2oo2 
1020cii2oi  -  4356c  11021  —  1 16c  10301  +  672c  102 11  +  4032cio2O2  +  672cioi2i  —  1228c  10031  + 


4032c  10022-^  35532c  10004  ~392co40oi  116c, .3101  -  1228cf,3oi  1 -^4704co3oo2- 174co220i  + 
672co21  1 1  +4032co21o2  ~  1842co202i  4032co2oi2  ~  21  ]68co2003  “  1  lOcoisoi  ■^072coi2i  1 

4032coi202  +  672coi  121  “  l228coio3i  4032cuio22 35532coioo4  27coo40i  +  448coo3i  1  “ 
324coo302  "♦'672coo22i  +  4032coo212  16542crio2Ci3  ^  448cooi31  4032cooi22  “  251cooo4i  + 
3012cooo32  ~  8478cooo23)/ 70056 

110:  cnooi4  =  (- Hc^oooi -64c3iomi  - 64c3f, 11,1  ’  I6C3001 1 168c3ooo2-96c220oj -216c2i  101  + 
24c2ion  ^  144c2  1002“  90c 20201  24f  21,1 1 1  ^  1 44c 201 02 24c 2002 1  144c2oo  12“ 756c 20003“ 
64  c  13001  “ 216c  12101 24c i2Cii  1  “  144c  12002  ~  216c  1 1201  24c  1 1021  “64c  10301  ■+‘24c  1021 1  ■*" 
144cio202  24cioi2i  “  16cioo3i  “  144ri,xj22  +  1269ciooo4  ~  14co4ooi  ~  64co3ioi  + 

16co3011  168co3002  “  96co22Ci1  24C02111  “  144co2102  24Co202I  144co2012  “ 

756Co2003  “  64Coi301  “  24Coi211  f  144coi202  “  24Coil21  +  IOCqiOSI  +  144co1022  + 

1269coio04  “  14coo401  ■*"  16C00311  “  168coo302  24C0022I  “  144Co0212  “  756Co0203  + 

16cooi3i  +  144ccoi22  f  1269cooio4  +  6cooo4i  “  72cooo32  1044cooo23)/3771 
111:  Coooos  =  (ciooo4  ~  C01004  +  C00104  C00014)  ^4 

112:  C01112  =  (“9c4oooi  “24c3iooi  “  24c3oioi  “  24c3oon  IO8C30002  “  36c2200i  “  36c2iioi  “ 
36C21011  +  54C21002  “  36c 20201  “  36c 20111  “  54C20102  “  36c 20021  54C20012  “  486C20003  “ 
24ci3ooi  “  36ci2ioi  -  36ci2oii  54ci2002  ~  36cn2oi  “  36ciio21  “  24cio3oi  “  36cio2ii  + 

54Cio202  “  36ci0121  “  24Ci0031  +  54Ci0022  “  972ci0004  f  C0400l  +  16Co3101  +  16Co3011  “ 

12co3002  +  24co2201  “  84Co21  11  “  36co2102  +  24Co2021  “  36co2012  +  54Co2003  +  16Coi301  + 
64Coi211  “  36co1202  f  84Coi121  t"  16coi031  “  36coi022  ~  108coio04  +  C00401  +  16co0311  ~ 

12Co0302  24Co0221  “  36co0212  +  54Coo203  +  16coo131  “  36coo122  “  108Co01O4  +  C00041  “ 

12cooo32  +  54cooo23  “  108coooi4  ~  162coooos)/180 
113:  C10112  =  (-2C40001  -  24c3iooi  +  4c3oioi  +  4c3oon  +  24c3ooo2  “  36c2200J  ~  36c2iioi  ~ 
36C21011  +  54C21002  +  6C2020I  +  48C20in  “  9C20102  +  6C2002I  “  9C20012  “  108C20003  “ 

24ci3ooi  ~  36ci2ioi  “  36ci2oii  +  54ci2002  -  36cii2oi  ~  36cno2i  +  4cio3oi  +  48cio2ii  “ 

9c  10202  48cioi21  +  4Ci0031  “  9cio022  +  216cj0004  “  6Co4001  “  12C03101  “  12co3011  + 
72Co3002  “  18co2201  27Co2102  “  18Co2021  +  27co2012  “  324co2003  “  12co1301  +  27coi202  “ 

54CoiU2  “  12coi031  +  27coI022  +  648coio04  +  C0040l  +  16co03I  I  “  12Co0302  +  24Co0221  ~ 

36Co0212  +  54Co0203  +  16co0131  “  36Coo122  “  108coO1O4  +  C00041  “  12Co0032  +  54Coo023  “ 

108coooi4  “  162coooo5)/126 

114:  c  11012  —  (“2C40001  “  8^31001  —  12C30101  +  4C30011  +  24C30002  “  12C22001  “  36C2n01  + 
12C2i011  +  I8C21002  “  I8C202OI  +  27C20102  +  6C2002I  “  9C200I2  “  108C20003  “  8Ci3001  ~ 
36ci2ioi  +  12ci2oii  I8C12002  “  36cji2oi  +  12ciio2i  -  12cio30i  +  27cio202  “  54cion2  + 
4Ci0031  “  9Cio022  +  216ci0004  “  2Co400I  “  12co3101  +  4Co301]  +  24C03002  ~  18C02201  + 
27Co2I02  +6co202I  -  9C02012  “  108co2OO3  ~  12Coi301  +  27co1202  “  54coin2  +  4Coi031  “ 
9Coi022  +  216Coi004  “  3Co0401  +  36Co0302  “  162Coo203  '*■  324Coo104  +  C0004I  “  12Co0032  + 
54cooo23  “  108coooi4  “  162coooo5)/72 

115:  C11102  =  (“C40001  “4C3ioo1  “  4C30I01  +  I2C30002  “6C2200I  -  12C21  101 -t-9C2i002“6C20201  + 
9C20102“54C20003- 4c  13001  -  1 2c  12101+ 9c  12002“  12c  11  201  -27c  1101 2  “4c  10301  +9^10202“ 
27Cioii2  +  108Cio004“C0400I  -4Co3101  +  12Co3002“ 600220 1  +9co2I02'54Co2003-4Co1301  + 
9coi202  “  27coiii2  +  108coioo4  “  C00401  12coo3o2  “  54coo203  +  108cooio4  ~  81coooo5)/9 
116:  Cl  1003  =  (C21002  +  C 12002  +  <^11102  +  Cl  1012) /'I 
117:  C10103  =  (C20102  +  C11102  +  C10202  +  cioii2)'^4 


118:  couo3  =  (ciuo2  +  C02102  cni202  coni2)/4 

119;  C10013  =  (C20012  +  C11012  ■+  cion2  cioo22)/4 

120:  C01013  =  (cuoi2  +  C02012  +  C0U12  +  coio22)/4 

121:  Coon3  =  (cioii2  +  coin2  C00212  corii22)/4 
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